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ABSTRACT 

Wc investigate the orbital evolution of a system of N mutually interacting stars on 
initially circular orbits around the dominating central mass. Wc include perturbative 
influence of a distant axisymmetric source and an extended spherical potential. In 
particular, we focus on the case when the secular evolution of orbital eccentricities is 
suppressed by the spherical perturbation. By means of standard perturbation methods, 
we derive semi-analytic formulae for the evolution of normal vectors of the individ- 
ual orbits. We find its two qualitatively different modes. Either the orbits interact 
strongly and, under such circumstances, they become dynamically coupled, precessing 
synchronously in the potential of the axisymmetric perturbation. Or, if their mutual 
interaction is weaker, the orbits precess independently, interchanging periodically their 
angular momentum, which leads to oscillations of inclinations. We argue that these 
processes may have been fundamental for the evolution of the disc of young stars 
orbiting the supermassive black hole in the centre of the Milky Way. 

Key words: methods: analytical - celestial mechanics - stars: kinematics and dy- 
namics - Galaxy: nucleus. 



1 INTRODUCTION 

Problem of dynamics in the perturbed Keplerian potential 
has been studied extensively throughout the history of celes- 
tial mechanics. Due to high attainable accuracy of observa- 
tional data, its primary field of application has always been 
the Solar System, which naturally infiuenced the selection 
of included perturbations. Ones of those widely considered 
are, due to their resemblance with the averaged motion of 
planets, axisymmetric gravitational potentials. 

The above problem has, however, also been investigated 
for systems with larger length scales, such as dense star clus- 
ters. In that case, the source of the Keplerian potential is of- 
ten represented by a supermassive black hole (SMBH) which 
is widely assumed to reside in the centres of such clusters. 
Axisymmetric perturbation is then either due to a secondary 
massive black hole (e.g. Ivanov, Polnarev & Saha 2005) or 
a gaseous disc or torus (e.g. Karas & Subr 2007). It turns 
out that in these systems, the secular evolution of individual 
stellar orbits is, beside the axisymmetric perturbation, also 
affected by a possible additional spherical potential. Such 
a potential may be generated by a stellar cusp or it can 
represent a post-Newtonian correction to the gravity of the 
central black hole. 
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In this paper, we extend the analyses of previous au- 
thors by means of standard tools of celestial mechanics. Our 
main aim is to incorporate mutual interaction of stars on 
nearly-circular orbits around the dominating central mass 
whose potential is perturbed by a distant axisymmetric 
source and an extended spherical potential. We apply our 
resul t s to the observed system of young stars (IGenzel et al 



resul t s to the observed system ot young stars 
20031: iGhez et all l2005l : iPaumard et"alll2006l : 



Bartko et al 



bood. (201011 orbiti n g the SMBH of rnass M . 4 x 10*^ Mr 
dGhez et all l2003l: lETsenhauer et all l2005l : iGillessen et~aLl 
l2009al lbl: lYelda et al.|]20ld ') in the "centre of the Milky Way 
As an axisymmetric perturbation to its gravity we consider 
a massive molecular torus (the so-called circumnuclear disc; 
CND) which is located at rad ius i?cND ~ 1-8 pc from the 
centre (|Christopher et al]|2005l ) . Finally, we co nsider gravity 
of a r oughly spherical cu s p of late-type s tars (|Genzel et al.l 
I2OO3I : ISchodel et al.|[2007l : Ido et al.H2009l ) which is believed 
to be present in this region, as well. Within this context, 
we broaden the analysis of our previous paper (Haas, Subr 
& Kroupa 2011) where we have studied the dynamical evo- 
lution of this kind of system purely by means of numerical 
AT-body calculations. In particular, we now develop a simple 
semi-analytic model which naturally explains key features of 
our prior results. 

The paper is organized as follows. In the theoretical 
Section [2I we first discuss the influence of the spherical per- 
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turbative potential upon the stellar orbits (Section l2.1|l . This 
allows us to separate the evolution of eccentricity from the 
rest of the problem and, subsequently, to formulate equa- 
tions for the evolution of inclinations and nodal longitudes 
(Section 12. 2|) . In Section [S] we present an example of the 
orbital evolution of a stellar disc motivated by the configu- 
ration that is observed in the Galactic Centre. We conclude 
our results in Section |4] 



2 THEORY 

To set the stage, we first develop a secular theory of orbital 
evolution for two (later in the section generalized to multi- 
ple) stars orbiting a massive centre, the SMBH, taking into 
account their mutual gravitational interaction and perturba- 
tions from the spherical stellar cusp and the axisymmetric 
CND. The CND is considered stationary and its model is 
further simplified and taken equivalent to a ring at a certain 
distance from the centre. It should be, however, pointed out 
that generalization to a more realistic structure, such as thin 
or thick disc, is straightforward in our setting but we believe 
at this stage it would just involve algebraic complexity with- 
out bringing any new quality to the model. In the same way, 
the stellar cusp is reduced to an equilibrium spherical model 
without involving generalizations beyond that level. For in- 
stance, an axisymmetric component of the stellar cusp may 
be effectively accounted for by the CND effects in the first 
approximation. 

We are going to use standard tools of classical celes- 
tial mechanics, based on the first-order secular solution us- 
ing the perturbation methods (see, e.g., Morbidelli 2002 or 
Bertotti, Farinella & Vokrouhlicky 2003 for general discus- 
sion). In particular, the stellar orbits are described using a 
conventional set of Kepler's elements which are assumed to 
change according to Lagrange equations. Since we are in- 
terrested in a long-term dynamical evolution of the stellar 
orbits we replace the perturbing potential (or potential en- 
ergy) with its average value over one revolution of the stars 
about the centre, which is the proper sense of addressing our 
approach as secular. In doing so, we assume there is no or- 
bital mean motion resonance between the two (or multiple) 
stars. As an implication of our approach, the orbital semi- 
major axes of the stellar orbits are constant and information 
about the position of the stars in orbit is irrelevant. The sec- 
ular system thus consists of description how the remaining 
four orbital elements, eccentricity, inclination, longitude of 
node and argument of pericentre, evolve in time. This is 
still a very complicated problem in principle, and we shall 
adopt simplifying assumption that will allow us to treat the 
eccentricities and pericentres separately (Section I2.1|l and 
leave us finally with the problem of dynamical evolution 
of inclinations and nodes (Section 12. 2p . Note this is where 
our approach diverges from typical applications in planetary 
systems, in which this separation is often impossible. 



describe under which conditions we may assume they stay 
small to the point we could neglect them. Note this is not an 
obvious conclusion because axially symmetric systems (such 
as a perturbing massive ring) have been extensively studied 
in planetary applications and it has been shown that non- 
conservation of the total orbital angular momentum may 
lead to a large, correlated variations of eccentricity and in- 
clination even if the initial eccentricity is arbitrarily small. 
This is often called Kozai secular resonance as a tribute to 
a pioneering work of Kozai (1962) (see also Lidov 1962). In 
what follows we describe conditions under which this process 
is inhibited in our model. 



2.1.1 Stellar cusp potential 

We start with our assumption about the potential energy 
of a star of mass m in the spherical cusp of the late-type 
stars surrounding the centre. Considering a general power- 
law radial density profile of the cusp, p (r) oc r~" , we have 
the potential energy 



7^c = - 



GmMc 



13 Rem \Rc 



(1) 



where /3 — 2 — a, the cusp mass within a scale distance -Rcnd 
is denoted Mc and G stands for the gravitational constant. 
According to the averaging technique, we shall integrate the 
potential energy ([!]) over one revolution about the centre 
with respect to the mean anomaly I, 



dl 7^c 



which yields 

1 GmA/c 



7^, 



2tv PRc^ 



dl 



r\/3 
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a/ 



(2) 



(3) 



where a and e are semi-major axis and eccentricity of the 
stellar orbit, r = a (1 — ecosu) and u — esinii — I. After an 
easy algebra, we obtain 

J{e,l3), (4) 



— _ GmMc 

/Vc 



/3i?CND V^CND 



where 



J(e,/3) = - / du (1 



with the coefficients obtained by recurrence 



a„ 



1 - 



3-f /3 



2(n-fl) 



1 - 



2 + li 



2(71+1) 



(5) 



(6) 



and an initial value ai = /3 (1 + /3) /4. For the purpose of 
our study, we further set /3 = 1/4 which corresponds to the 
equilibrium model worked out bv lBahcall fc Wol3 (|l976l '). 



2.1 Confinement of eccentricity 

In this section we discuss our assumptions about eccentricity 
and pericentre evolution. For this moment, we drop the mu- 
tual interaction of stars from our consideration. We assume 
that the initial stellar orbits have small eccentricity and we 



2.1.2 Circumnuclear disc/ring potential 

In the case of perturbation of orbits well below the radius of 
the CND, we limit ourselves to account for the quadrupole- 
tide formulation (e.g. Kozai 1962; Morbidelh 2002). Oc- 
tupole or higher-multipole corrections are possible (e.g. in 
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Figure 1. Isolincs of the conserved potential function 1^ = 
from equation (|9]l for two different values of the mass ratio 
11 = Mc/McND- 0.01 at the top panel, 0.1 at the bottom panel. 
The Kozai integral value is c = cos(70°), corresponding to 70° 
inclination circular orbit. The orbit has been given semi-major 
axis a = 0.06 -Rcnd for sake of definiteness. The origin e = is a 
stationary point of the problem but in the first case it is unsta- 
ble, while in the second case it becomes stable. The thick isoline 
in the top panel is a separatrix between two different regimes of 
eccentricity and pericentre evolution. 



fact Kozai himself gives explicit terms up to degree 4; see 
also Yokoyama et al. 2003) but they do not change the con- 
clusions as long as the parameter a/ Rcnd is small enough. 
This is the regime that interests us most. 

Given the axial symmetry of the mass distribution of 
the perturbing ring, the resulting averaged interaction po- 
tential energy of a particle in the tidal field of the CND (see 
Kozai 1962) 



7?.CND = — 



GmMcND / a 
16J?CND V-^CN 



-I- 15e^ sin^ I cos 2lj 



-] [(2-f 3e^) (Scos^J- 1) 



(7) 



does not depend on longitude of node but depends on 
other orbital elements of the stellar orbit - eccentricity e, 
inclination I and argument of pericentre lu. 

As a direct consequence, c = V 1 — cos / is the first 
('Kozai') integral of motion, which conveniently allows to 
eliminate inclination dependence in 7?.cnd, depending then 
on the eccentricity and argument of pericentre only. Since 
7?.CND is a conserved quantity in the secular (orbit-averaged) 




0.5 



Figure 2. Individual lines show a critical inclination (ordinate) 
at which Kozai resonance onsets for a given value of mass ratio 
H = Mc/McND (abscissa) for different values of orbital semi- 
major axis a ranging from 0.03 -Rcnd (left) to 0.3-RcND (right) 
with the step of 0.03 -Rcnd- When /i = 0, the critical angle is 
~ 39.2° ('the Kozai limit') independently from a. 



problem, the isolines 7?.cnd ~ C provide insights in the fun- 
damental features of the dynamical evolution of both e and 
ij. This approach has been used by Kozai to discover two 
modes of topology of these isolines: (i) when c > \/3/5 the 
T^-CND = C isolines are simple ovals about origin which is the 
only fixed point of the problem but (ii) for c ^ ^3/5 they 
become more complicated with a separatrix curve emerging 
from the origin and two new fixed points exist at nonzero ec- 
centricity and pericentre argument values 90° and 270° . The 
latter case occurs whenever the initial inclination is larger 
than ~ 39.2°, sometimes called the Kozai limit. The impor- 
tant take-away message is that the circular orbit is no more a 
stable solution for high-inclination orbits in the model of ex- 
terior ring/disc perturbation. Initially circular orbits would 
be driven over a Kozai timescale 



Tk 



M. 



McND aVGM.a 



(8) 



to a very high eccentricity state. Unavoidable stellar scatter- 
ing processes would in a short time destabilize an initially 
coherent stream of objects near the centre. 



2.1.3 Combined perturbation 

We now consider combined effect of the stellar cusp and the 
CND potentials on the long-term orbital evolution of the 
stellar orbit. The total, orbit-averaged potential 



7^ = Tic + TlcND 



(9) 



still obeys axial symmetry, being independent on the nodal 
longitude. The picture, however, may be modified with re- 
spect to the case of solely ring-like perturbation. Considering 
the cusp of the late-type stars whose potential is approxi- 
mated with Q , we find that the two types of topologies of 
the TZ = C isolines persist (see Fig. [T| but the onset of the 
circular-orbit instability depends now on two parameters, 
namely c and fx = Mc/Mcnd- A nonzero mass of the stellar 
cusp stabilizes small eccentricity evolution and the critical 



4 J. Haas, L. Subr and D. Vokrouhlicky 



angle is pushed to larger values. For large enough n, the sta- 
bility of the circular orbit is guaranteed for arbitrary value 
of c and hence orbits of an arbitrary inclination with respect 
to the CND symmetry plane. This is because the effects of 
the stellar cusp potential make the argument of pericen- 
tre circulate fast enough (significantly faster than the Kozai 
timescale), preventing thus secular increase of the eccentric- 
ity. An initially near-circular orbit maintains a very small 
value of e showing only small-amplitude oscillations. Fig. [2] 
shows critical inclination values, for which the circular orbit 
becomes necessarily unstable as a function of fi and a/_RcND 
parameters (note the later factorizes out from the analysis 
when /i = 0). Importantly, there is a correlation between 
H and a/RcND below which circular orbits of an arbitrary 
inclination are stable; for instance, data in Fig. [2] indicate 
that for /i = 0.1 any circular orbit with a < 0.12_Rcnd is 
stable. 

In conclusion, we observe that having enough mass in 
the late-type stellar cusp may produce strong enough per- 
turbation to maintain small eccentricity of an initially near- 
circular orbit. With that said, we find it reasonable to make 
an important simplification within our analytic approach to 
the system of two (multiple) stars. Namely, we will further 
consider the stellar orbits to be circular during the whole 
evolution of the system. This prevents (together with the as- 
sumption of well separated orbits with constant semi-major 
axes) close encounters of the stars. In this case only, and 
under the assumption that there are no orbital resonances 
among the individual stars, the mutual interaction of the 
stars may be reasonably considered as a perturbation to the 
dominating potential of the SMBH. As we demonstrate in 
the next sections, this simple treatment provides useful in- 
sights into the evolution of the young-stream orbits even if 
they are generally non-circular. 



2.2 Orbital evolution of circular orbits 

Having discussed our assumptions about semi-major axes, 
eccentricity and pericentre of the stellar orbits, we may now 
turn to description of the evolution of the two remaining 
orbital elements - inclination and nodal longitude. We start 
with a model of two interacting stars and later generalize 
it to the case of an arbitrary number of stars. The major 
leap-forward in the model is that we now take into account 
also mutual gravitational effects of the two stars. On the 
contrary, note that the orbit-averaged potential energy Q 
of the late-type stellar cusp depends on the semi-major axis 
and eccentricity only, and thus does not influence evolution 
of inclination and node. For that reason it drops from our 
analysis in this section. 

The interaction potential energy TZi{r, r') for two point 
sources of masses m and m' at relative positions r and r' 
with respect to the centre readfl 



7^i(r,r') = 



Gmm' 



^ a Pi (cos S) 



(10) 



^ Note that equation IjlOj l provides the interaction energy as it 
appears in the equation of relative motion of stars with respect 
to the centre. Henceforth, the perturbation series start with a 
quadrupole term [1 = 2). 



where Pi{x) are Legendre polynomials, cos S = r-r' /rr' and 
a = r' jr. The series in the right-hand side of equation (|10|l 
converge for r' < r. Since we are going to apply (|10|l to 
the simplified case of two circular orbits, we may replace 
distances r and r' with the corresponding values of semi- 
major axis a and a', such that a = a' /a now (note that the 
orbit whose parameters are denoted with a prime is thus 
assumed interior). The averaging of the interaction energy 
over the uniform orbital motion of the stars about the centre, 
implying periodic variation of S, is readily performed by 
using the addition theorem for spherical harmonics. This 
allows us to decouple unit direction vectors in the argument 
of the Legendre polynomial Pe and easily obtain the required 
average of TZi over the orbital periods of the two stars. After 
a simple algebra we obtain 



7^i = * (a,n • n ) , 



(11) 



where n — [sin / sin f2, — sin / cos f2, cos 7]^ and n' — 
[sin 7' sin f2', — sin J' cos cos J']^ are unit vectors normal 
to the mean orbital planes of the two stars, and 



*(C,^) = E[ft(0)]'c'ft W 



(12) 



As expected, the potential energy is only a function of: (i) 
the orbital semi-major axes through dependence on a and 
a, and (ii) the relative configuration of the two orbits in 
space given by the scalar product n ■ n' . Note also that the 
series in (|12|l contain only even multipoles I {Pe.{Q) = for 
£ odd) and that they converge when C < 1- However, a spe- 
cial care is needed when ^ is very close to unity, thus the 
two stellar orbits are close to each other, when hundreds 
to thousands terms are needed to achieve sufficient accu- 
racy. Still, we found it is very easy to set up an efficient 
computer algorithm, using recurrent relations between the 
Legendre polynomials, which is able to evaluate (|12p and its 
derivatives. In practice, we select a required accuracy and 
the code truncates the series by estimating the remained 
terms. In fact, since our approach neglects small eccentric- 
ity oscillations of the orbits we are anyway not allowed to 
set C, = a = a' /a arbitrarily close to unity. Theoretically, we 
should require 



a < 1 



m + m 
3A/. 



1/3 



(13) 



by not letting the stars approach closer than the Hill radius 
of their mutual interaction. In the numerical examples we 
present below, this sets an upper limit a < 0.98. 

The formulation given above immediately provides po- 
tential energy of the star-CND interaction. In this case the 
stellar orbits are always interior to the CND with symmetry 
axis suitably chosen as the unity vector in the direction 
of the z-axis of our reference system. Unlike in Section [2. 1.21 
we restrict now to the case of circular orbit of the star but 
at the low computer-time expense we may include all multi- 
pole terms till specified accuracy is achieved. As a result the 
orbit-averaged interaction energy with the exterior stellar 
orbit is given by 



GmMcND 

7?CND 



* (a/iicND,cos7) , 



(14) 
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and similarly for the interior stellar orbit: 
Gm'McND 



7^, 



* (a7-RcND,cos/') 



(15) 



The total orbit-averaged potential energy perturbing motion 
of the two stars is then given by superposition of the three 
terms: 



7^ = 7^i + 7^cND + iZc 



(16) 



Recalling that semi- major axis values are constant, eccen- 
tricity set to zero and thus argument of pericentre undefined, 
we are left to study dynamics of inclination I and I' and lon- 
gitude of node Q and il' values. Lagrange equations provide 
(see, e.g., Bertotti et al. 2003) 



mna^ 9 cos/ ' 

1 dn 



(17) 
(18) 



dcosJ _ 1 dTZ dQ 1 dTZ 

dt mna^ dO, ' dt 

dcosJ' _ 1 dn dO.' 

dt m'n'a''^ do,' ' dt m'n'a''^ d cos I' 

where n and n denote mean motion frequencies of the two 
stars. Note the particularly simple, quasi-Hamiltonian form 
of equations (|17p and (|18[l . They can also be rewritten in a 
more compact way using the normal vectors n and n' to the 
respective orbit, namely 



dn d 
-T7 = nx — 
dt on 



n 



dn' 



d 



= n X „ 
dt dn' 



7^ 



Inserting here TZ from (|16|l . we finally obtain 

dn , , , 

— = oji[n xn) + CJCND (n x e^) , 

dn' 



dt 
where 

OJl = - 
UJi^ - 



LUi in X n) + CJCND yn x e.^ 



m 

m7 

m 
IT. 

, / A/cND \ 



a, n ■ n') , 

(a/i?cND,nz) , 
•^x [a / RcND , n'^) 



(19) 
(20) 

(21) 
(22) 

(23) 
(24) 
(25) 
(26) 



Note the frequencies in (|23|) to (|26p depend on both n and 
n through their presence in the argument of 

d 



(27) 



which breaks the apparent simplicity of the system of equa- 
tions (dll and ((22)1 . 

The coupled set of equations (|21[l and (|22p acquires sim- 
ple solutions in two limiting cases. First, when m = m' — 
(i.e. mutual interaction of stars is neglected) the two equa- 
tions decouple and describe simple precession of n and 
n' about e, axis of the inertial frame with frequencies 
— ojcND cos / and — oj^nd cos /'. The sign minus of these fre- 
quencies recalls that the orbits precess in a retrograde sense 
when inclinations are less than 90° and vice versa. Both in- 
clinations / and /' are constant. In the second limit, when 




180° 



Figure 3. Isolines of the 7?. = C2 integral in the / or /' vs. 
Af2 space. For sake of example we use orbits of two equal-mass 
stars (m' = m) with semi-major axes a' = 0.04i?cND and a 

■ The mass of the CND is set to AfcND = 0.3 Af.. 

The individual lines correspond to different values of stellar mass: 
m = 5 X 10"'^ M, (curves 1), m = 2 X 10"'' M, (curves 2), m = 
5 X lO""" M. (curves 3), and m = 9 X 10"^ M, (curves 4). Both 
orbits have been given 70° inclination at Af2 = 0° (i.e. initially 
coplanar and inclined orbits). Solid lines show inclination /' of the 
inner orbit, 'the mirror-imaged' dashed lines describe inclination 
I of the outer orbit. 



AfcND = (i.e. the circumnuclear torus is removed) the 
equations (|2ip and (I22p obey a general integral of total an- 
gular momentum conservation 

mn + m'a^^^ n — K . (28) 

Both vectors n and n' then precess about K with the same 
frequency 



m + m'a ' {n ■ n') 



m'a^^^ ^m? -f m'^a -f 2mm' a^l'^ (n • n') 



(29) 



keeping the same mutual configuration. In particular, ini- 
tially coplanar orbits (i.e. n and n' parallel) would not 
evolve, which is in agreement with intuition. 

Unfortunately, we were not able to find analytical solu- 
tion to the (|2ip and (|22p system except for these two situ- 
ations described above. Obviously, it can be always solved 
using numerical methods as we shall discuss in Section r2.2.2l 



2,. 2.1 Integrals of motion 



In general, equations (|21|l and (|22p have only two first inte- 
grals. Our assumptions about the circumnuclear torus mass 
distribution still provide a symmetry vector e^. Thus, while 
the total angular momentum K is no more conserved now, 
its projection onto is still an integral of motion 

m cos I + m'a^^^ cos l' — Ci — Kz ■ (30) 

Because m, m' and a are constant, equation (|30p provides a 
direct constraint of how the two inclinations I and /' evolve. 
In particular, one can be expressed as a function of the other. 

The quasi-Hamiltonian form of equations (|17p and ([18} 
readily results in a second integral of motion 

7?. (cos /, cos /', n • n') = C2 . (31) 
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Figure 4. Evolution of the system of two stars in the eompound potential of the eentral SMBH, spherical stellar eusp and axisymmetric 
CND. Solid lines represent solution of two-body equations II2H and ([22}, while the dashed lines show result of the direct numerical 
integration of the equations of motion. In each panel, upper and lower lines correspond to the inner and outer star, respectively. Common 
parameters for both examples are the same as in Fig. \3\ in the upper panels, we set m = m' = 9 X 10~® Mm , while in the lower ones 
m = m' = 5 X 10-6 M. . 



The list of arguments in TZ, as explicitly provided above, 
reminds that it actually depends on: (i) the inclination val- 
ues / and /', and (ii) the difference AQ — — il' of the 
nodal longitudes of the two interacting orbits. Using (|30p . 
the conservation of TZ thus provides a constraint between 
the evolution of I and ASl (say). While not giving a solu- 
tion of the problem, the constraint due to combination of 
first integrals (|3Up and (|3ip can still provide useful insights. 

Fig.[3]illustrates how the first integrals help understand- 
ing several features of the orbital evolution for two interact- 
ing stars at distances a' — 0.04i?cND and a = 0.05 7?cnd. 
For sake of simplicity we also assume their mass is equal, 
hence m' = m, and the CND has been given mass AfcND = 
0.3 M,. Data in this figure show constrained evolution of 
orbital inclinations /' (solid lines) and / (dashed lines) as 
a function of nodal difference AJl. The two orbits were as- 
sumed to be initially coplanar (An = 0°) with an inclination 
oi I' = I — 70° . A set of curves correspond to different val- 
ues of stellar masses, from small (1) to larger values (4), 
which basically means increasing strength of their mutual 
gravitational interaction. 

First, conservation of the Cz-projected orbital angular 
momentum, as given by equation (|30[) . requires that increase 
in /' is compensated by decrease of /. This results in a near- 



mirror-imaged evolution of the two inclinations. Using the 
first equation of (|17|) . one finds 

S-ii^7S^^^'^("-"')*^("'--')- (32) 

which straightforwardly implies that the outer stellar orbit is 
initially torqued to decrease its inclination while the inner 
orbit increases its inclination. This is because initially n ■ 
n' ~ 1, and is positive, and at the same time, 

precession of the nodes is dominated by interaction with the 
CND which makes the outward orbit node to drift faster 
(and hence f2 — f2' is negative). 

Second, Fig. [3] indicates there is important change in 
topology of the isolines TZ — C2 as the stellar masses over- 
pass some critical value (about 8.5 x 10~^ Mm in our ex- 
ample). For low-mass stars their mutual gravitational in- 
teraction is weak letting the effects of the CND dominate 
(curve 1). The orbits regularly precess with different fre- 
quency, given their different distance from the centre, and 
thus AQ acquires all values between —180° and 180°. The 
mutual stellar interaction produces only small inclination os- 
cillation. As the stellar masses increase (curves 2 and 3) the 
inclination perturbation becomes larger. For super-critical 
values of m (curve 4) the isolines of constant TZ become 
only small loops about the origin. This means that ASl is 
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Figure 5. Evolution of tiie system of four stars in the compound potential of the eentral SMBH, spherical stellar cusp and axisymmetric 
CND. The stellar orbits form two couples. In both of them, the orbits have similar semi-major axes in order to mimic the system shown 
in Fig.|4] In each panel, upper and lower lines correspond to the inner and outer couple, respectively. The individual semi-major axes are 
for both examples set to ai = 0.0373 RcnDj 12 = 0.0408 RcnDj 13 = 0.0478 RcnDi 14 = 0.0511 Rcnd- The other common parameters 
for both examples are the same as in Fig.fS] in the upper panels, we set mi = m2 = = 7714 = 4.5 X 10~^ Mm, while in the lower ones 
mi = m2 = m:i = 7714 = 2.5 X lO^'' A/.. 



bound to oscillate in a small interval near origin and inclina- 
tion perturbation becomes strongly damped. Put in words, 
the gravitational coupling between the stars became strong 
enough to tightly couple the two orbits together. Note that 
they still collectively precess in space due to the influence of 
the CND. 



2.2.2 Numerical solutions 

In order to solve equations (|21|l and (|22|) numerically, we 
adopt a simple adaptive step-size 4.5th-order Runge-Kutta 
algorithm. Let us mention that our implementation of this 
algorithm conserves the value of both integrals of motion Ci 
and C2 with relative accuracy better than lO"*". 

Two sample solutions are shown in Fig. 2] The upper 
panels represent evolution of two orbits with coupled pre- 
cession which corresponds to the curve 4 in Fig. [S] while in 
the bottom panels we consider the case of lower-mass stars, 
whose orbits precess independently. This later mode corre- 
sponds to the curve 3 in Fig. [S] Beside the solution of the 
equations for mean orbital elements, we also show results 
of a full-fledged numerical integration of the particular con- 
figuration in the space of classical positions and momenta 



{r,r';p,p'). Both solutions are nearly identical, which con- 
firms validity of the secular perturbation theory used in this 
paper. 

For sake of further discussion we find it useful to com- 
ment in a little more detail on the case of two, nearly in- 
dependently precessing orbits (bottom panels on Fig. [Jjl . In 
this case, the precession frequencies of the outer and inner 
star orbits are given by ojcnd and i^cnd in equations p5|l 
and (|26|l . When truncated to the quadrupole (£ = 2) level, 
sufficient for the small value of a/RcND, one has for the 
outer star orbit 

dfi 3 cos J 
d^-~4^' 

where Tk is given by ((8]). A similar formula holds for the 
inner star orbit denoted with primed variables. As seen in 
Fig. O and understood from the analysis of integrals of mo- 
tion in Section r2.2.1l period of the evolution of the system of 
the two orbits is given implicitely by the difference of their 
precession rate: fi(rchar) — i^'(?char) = 27r. The secular rate 
of nodal precession in (|33|) is not constant because the mu- 
tual gravitational interaction of the stars makes their orbital 
inclinations oscillate. However, in the zero approximation we 
may replace them with their initial values, I = I' = Iq which 
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Figure 6. Evolution of tiie initially thin stellar disc of 100 stars in the compound potential of the central SMBH, spherical stellar cusp 
and axisymmetric CND. The values of orbital semi-major axes in the disc range from 0.02i?cND to 0.2RcND ^^nd their distribution 
obeys dN <x a^^da. The stellar masses are all equal with m = 5 X 10~^ M, while the mass of the CND is set to AifcND = 0.3 M,. Initial 
inclination Iq of all the orbits with respect to the CND equals 70° . 



gives an order of magnitude estimate 



Stt 



3 cosJo 



1 



1 



(34) 



For the solution shown in the lower panels of Fig. |4j for- 
mula (|34p gives Tchar ~ 460 Myr, in a reasonable agreement 
with the observed period of « 140 Myr. When the orbital 
evolution is known (being integrated numerically), more ac- 
curate estimate can be obtained considering mean values of 
the inclinations 



T ~ 

char — ^ 



COS / cos I 



Tk 



T' 



(35) 



For the case of the solution of the lower panel of Fig. |4l with 
7 ^ 60° and 7' a; 80°, formula ^ gives Tchar ~ 120 Myr. 



2.2.3 Generalization for N interacting stars 

The previous formulation straightforwardly generalizes to 
the case of A'^ stars orbiting the centre on circular orbits 
with semi-major axes at (fc = 1, . . . , A'^). This is because the 
potential energies of all pairwise interactions built the total 



TCi = -- — w (aki,nk. ■ ni) 



(36) 



where is the mass of the k-th star, a^i = min(afe,ai), 
aki = min(afe, ai)/max(at;, a;) and rife is the normal vector 
to the orbital plane of the fc-th star. Similarly, interaction 
with the CND is simply given by 

Hem = - Gmfc A/cND ^ (afc/_RcND, nfc ■ e^) . (37) 

k 

The total potential energy of perturbing interactions is 
TZ = TZi+ TZcNu , (38) 
and the equations of orbital evolution now read 



drik d 

—r— = rife X — , 

dt duk \mknka% 



(39) 



for k = 1, . . . , A'^ (rife is the frequency of the unperturbed 
mean motion of the fc-th star about the centre). Their first 
integrals then can be written as 



^mfenfeOfe (jife ■ ez) = Ci = 



and 



7^ = Ca 



(40) 



(41) 



Due to mutual interaction of multiple stars, solutions 
of equations H39p represent, in general, an intricate orbital 
evolution, whose course is hardly predictable as it strongly 
depends upon the initial setup. Our numerical experiments 
show, however, that it is still possible to identify several 
qualitative features which remain widely valid. For instance, 
a group of orbits with small separations may orbitaly couple 
together and effectively act as a single orbit in interaction 
with the rest of the stellar system. 

This is demonstrated in Fig. [S] which shows two sample 
solutions of equations (|39|l for a system of two such groups. 
For sake of clarity, each group consists only of two orbits. 
Individual semi-major axes are, for both solutions, set to 
ai = 0.0373 i?cND, a2 = 0.0408 TJcnd, as = 0.0478 i?cND, 
04 = 0.0511 -RcND in order to mimic the two-orbits models 
from Fig. |4] For the same reason, all the individual masses 
are considered equal, mi = m2 = ms = m^, and set to 
2.5 X 10~® M, in the lower panels, while for the upper pan- 
els we assume 4.5 x 10~^ M. . The other parameters remain 
identical to the case of the two-orbits models. As we can 
see (cf. Figs [S] and 131, the dynamical impact of each coupled 
pair of orbits upon the rest of the stellar system is equiv- 
alent to the effect of the corresponding single orbit if both 
the total mass and semi-major axis of the pair are appropri- 
ate. The individual orbits within each pair then naturally 
oscillate about the single-orbit solution according to their 
mutual interaction. This conclusion remains valid even in 
more complicated systems as we shall show in the next sec- 
tion. 
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3 APPLICATION TO THE YOUNG STELLAR 
SYSTEM IN THE SGR A* REGION 

In order to illustrate the complexity of solutions of equa- 
tions (|39p . let us now analyze the evolution of a system which 
contains an initially thin stellar disc with a distribution of 
semi-major axes of the orbits dA^ oc a^^da. As we can see in 
Fig. m the oscillations of the orbital inclinations no longer 
have the simple patterns which we observed for the models 
discussed in the previous paragraphs. On the other hand, 
we still can identify a well defined group of orbits which co- 
herently change their orientation with respect to the CND. 
These orbits thus form a rather thin disc during the whole 
monitored period of time. It turns out that they represent 
the innermost parts of the initial disc where the separations 
of the neighbouring orbits are small enough for their mutual 
interaction to couple them together. 

The configuration considered in Fig. [6] roughly matches 
the main qualitative features of an astrophysical system 
which is observed in the centre of the Milky Way. It contains 
a group of early-type stars orbiting the SMBH on nearly Ke- 
plerian orbits. Observations suggest that about one half of 
them form a coherently rotating disc-like st ructure with es- 
timated surface den s ity profile S oc (|Paumard et al.l 
I2OO6I : ILu et~al]|2009l : iBartko et al.ll2009l ) which implies the 
above considered distribution of semi-major axes. The rest of 
the early-type stars then appear to be on randomly oriented 
orbits. Both the origin and observed configuration of these 
stars represent rather puzzling questions. Due to strong tidal 
field of the SMBH, it is impossible for a star to be formed 
in this region by any standard star formation mechanism. 
On the other hand, as the observed stars are assumed to 
be young, no usual transport mechanism is efficient enough 
to bring them from farther regions, where their formation 
would be less intricate, within their estimated lifetime. One 
of the most promising scenarios of their origin thus consid- 
ers formation in situ, via fragmentation of a self-gravitating 
gaseous disc ([Levin fc Beloborodwl 120031 ). However, since 
this process naturally forms stars in a single disc-like struc- 
ture, it does not explain the origin of the stars observed 
outside the disc. Hence, in order to justify the in-disc sce- 
nario of the formation of the early-type stars in the Galactic 
Centre, some mechanism that may have dragged some of 
them out from the parent stellar disc plane is needed. 

In our previous paper (|Haas et al.ll20lil ). we have dis- 
cussed a possibility that all the early-type stars had been 
born in a single disc which has been, subsequently, partially 
disrupted by the gravity of the CND. We have considered 
the same configuration of the sources of the gravitational 
field as in the current paper and followed the evolution of 
the disc by means of direct A*'-body integration. We have ob- 
served coherent evolution of the inner dense part of the disc 
which exhibited a tendency to increase its inclination with 
respect to the CND. On the other hand, most of the orbits of 
the outer parts of the initially coherently rotating disc pre- 
cessed independently due to the influence of the CND and, 
consequently, became detached from the parent structure. 
This behaviour is in accord with the analysis presented in 
the current paper. 

Furthermore, we can now calculate the order of magni- 
tude chara c teristi c time-scale for the 'canonical' model of 
iHaas et all (|201ll ) whose system parameters read: A/. = 



4xlO'^M0, iJcND = l.Spc, McND =0.3M., Mc =0.03M., 
and Jo = 70°. In order to determine the rough time estimate, 
we use formula (|34p . As this formula has been derived for 
a system of two stars, we replace the stellar disc with two 
characteristic particles at certain radii a', a in the sense of 
Section [2. 2. 3 1 For this purpose, let us divide the stars in the 
disc into two groups according to their initial distance from 
the centre and define a' and a as the radii of the orbits of the 
median stars in the inner and outer group, i.e. a' — 0.06 pc 
and a — 0.23 pc. Inserting these values into formula (|34p . we 
obtain Tchar ~ 37 Myr for the 'canonical' model. This value 
is in order of magnitude agreem ent with the estimated age of 
the early-type stars, ~ 6 Myr jPaumard et al.ll2006l ). since 
the core of the disc reaches its maximal inclination with re- 
spect to the CND already after a fraction of period Tchar as 
can be seen in Figs. [5] and [6] 

Let us emphasize that t he results reported in our previ- 
ous paper (|Haas et al.ll201ll ) have been acquired by means of 
full-fledged numerical integration of equations of motion. As 
a consequence, both the eccentricities and semi-major axes 
of the individual stellar orbits in the disc have been nat- 
urally undergoing a significant evolution due to two-body 
relaxation of the disc. Moreover, our prior numerical com- 
putations have also confirmed that results similar to those 
obtained for the 'canonical' model are valid for a wide set of 
models with different system parameters, including the case 
with zero mass, Mc, of the spherical cusp of the late- type 
stars. In the later case, the orbital eccentricities and incli- 
nations within the stellar disc are subject to high-amplitude 
Kozai oscillations. In conclusion, it appears that the inner 
part of the disc may evolve coherently for a certain period of 
time even when we cannot assume neither zero nor small ec- 
centricity of the stellar orbits. We, therefore, suggest that 
also some of the key qualitative predictions of the semi- 
analytic theory developed in the current paper under the 
simplifying assumption of circular orbits may be carefully 
applied to more general, non-circular systems. 

Finally, let us mention that, in addition to the core of 
the disc, less significant groups of orbits with coherent sec- 
ular evolution may exist even in the outer parts of the disc 
if their separations are small enough. Our semi-analytic ap- 
proach thus admits possible existence of secondary disc-like 
structures in the observed young stellar s ystem which has in- 
dccd been discussed by several authors jCenzel et al.ll2003l : 
[Paumard et al.ii200^ : iBartko et al.ll2009f ). 



4 CONCLUSIONS 

We have investigated the secular orbital evolution of a sys- 
tem of mutually interacting stars on nearly-circular or- 
bits around the dominating central mass, considering the 
perturbative gravitational infiuence of a distant axisymmet- 
ric source and an extended spherical potential. Given the 
spherical potential is strong enough, we have shown that 
the secular evolution of initially circular orbits reduces to the 
evolution of inclinations and nodal longitudes. The spherical 
potential itself can then be factorized out from the outcom- 
ing momentum equations. Since we have not been able, in 
a general case, to solve the derived equations analytically, 
we have set up an integrator for their efficient numerical 
solution. The acquired results have then been, in order to 
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confirm their validity, compared witii tlie corresponding full- 
fledged numerical integrations in the space of classical posi- 
tions and momenta, showing a remarkable agreement. 

Some fundamental features of the possible solutions of 
the new equations can be understood by an analysis of the 
integrals of motion. In the case of the simplest possible sys- 
tem of two stars interacting in the considered perturbed po- 
tential, we have identified two qualitatively different modes 
of its secular evolution. If the interaction of the stars is weak 
(yet still non-zero), the secular evolution of their orbits is 
dominated by an independent nodal precession. Difference 
of the individual precession rates then determines the pe- 
riod of oscillations of the orbital inclinations. On the other 
hand, when the gravitational interaction of the stars is suf- 
ficiently strong (depending on their mass and the radii of 
their orbits), the secular evolution of the orbits becomes 
dynamically coupled and, consequently, they precess coher- 
ently around the symmetry axis of the gravitational poten- 
tial. Oscillations of the orbital inclinations are, in this case, 
considerably damped. 

We have further confirmed, by means of numerical in- 
tegration of the derived momentum equations, that the cou- 
pling of strongly interacting orbits is a generic process that 
may occur even in more complex A^'-body systems. In partic- 
ular, a subset of stars with strong mutual interaction evolves 
coherently and, as a result, its dynamical impact upon the 
rest of the A''-body system is similar to the effect of a single 
particle of suitable mass and orbital radius. 

As an example, we have investigated evolution of a disc- 
like structure that roughly models the young stellar system 
which is observed in the Galactic Centre. It has turned out 
that the semi-analytic work presented in this paper pro- 
vides a physical background for understanding of the pro- 
cesses discovered, by means of full A^'-body integration, in 
iHaas et al.l (|201ll V Namely, coupling of the strongly inter- 
acting stars from the inner parts of the disc leads to their 
coherent orbital evolution, which allows us to observe a disc- 
like structure even after several million years of dynamical 
evolution in the tidal field of the CND. Orientation of this 
surviving disc then inevitably changes towards higher incli- 
nation with respect to the CND, which is in accord with 
the observations. On the other hand, stellar orbits from the 
outer parts of the disc evolve individually, being gradually 
stripped out from the parent thin disc structure. Hence, it 
appears possible for the puzzle of the origin of the young 
stars in the Galactic Centre to be solved by the hypothesis 
of their formation via fragmentation of a single gaseous disc, 
as a lready suggested in Subr, Schovancova & Kroupa (2009) 
and lHaas etal] (|201ll ). 

Note that, beside the physical explanation of the pro- 
cesses observed in our previous work, the current approach 
would be, due to its low numerical demands, useful for ex- 
tensive scanning of the parameter space in order to confront 
our model with the observations more thoroughly. This is 
going to be a subject of our future work when more accu- 
rate observational data will be available. 

Finally, let us mention that our semi-analytic model has 
been developed under several simplifying assumptions. Most 
importantly, the torus CND has been considered stationary 
and the cusp of the late-type stars spherically symmetric. 
If any of these assumptions were violated, the results might 
be more or less affected. For example, a possible anisotropy 



of the cusp of the late-type stars due to chance alignment 
of some of its stars would break its spherical symmetry. In 
that case, the resulting gravitational torques might have a 
considerable imp act on the dynamica l evolu tion of the stellar 
disc as shown bv lKocsis fc Tremainel (|201ll ). However, since 
the current observational data do not show evidence for such 
violations, we may consider our model physically plausible. 
Moreover, the currently available data do suggest roughly 
perpendicular mutual orientation of the CND and the stellar 
disc, which is in accord with the predictions of both our 
numerical and semi-analytic model. We consider this as a 
supporting argument for our findigs. 
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